\documentclass[12pt]{ctexart}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage[UTF8]{ctex}
\usepackage{amsfonts,amssymb}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{booktabs}
\usepackage{listings}
\usepackage{geometry}
\usepackage[hidelinks]{hyperref}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{setspace}
\usepackage{float}   
\usepackage{svg}
\usepackage{url}
\usepackage{epstopdf}
\usepackage{tikz}

\geometry{a4paper}
\geometry{top=3cm} 
\geometry{bottom=3cm}
\geometry{left=2.5cm}
\geometry{right=2.5cm}
\definecolor{mygreen}{rgb}{0,0.6,0}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
\definecolor{mymauve}{rgb}{0.58,0,0.82}

\usepackage{listings}
\lstset{
 linewidth=1.1\textwidth,
 basicstyle=\small\ttfamily, % 设置字体族
  backgroundcolor=\color[RGB]{245,245,244},
 breaklines=true, % 自动换行
 keywordstyle=\bfseries\color{NavyBlue}, % 设置关键字为粗体，颜色为 NavyBlue
 morekeywords={}, % 设置更多的关键字，用逗号分隔
 emph={self}, % 指定强调词，如果有多个，用逗号隔开
    emphstyle=\bfseries\color{Rhodamine}, % 强调词样式设置
    commentstyle=\itshape\color{black!50!white}, % 设置注释样式，斜体，浅灰色
    stringstyle=\bfseries\color{PineGreen!90!black}, % 设置字符串样式
    columns=flexible,
     tabsize=4,
      escapeinside={\%*}{*)}, 
    numbers=left, % 显示行号在左边
    numbersep=2.5em, % 设置行号的具体位置
    numberstyle=\tiny, % 缩小行号
    frame=single, % 边框
    framesep=1em, % 设置代码与边框的距离
    extendedchars=false, %解决代码跨页时，章节标题，页眉等汉字不显示的问题  
    xleftmargin=2.5em,xrightmargin=2.5em, aboveskip=1em, %设置边距  
    columns=fullflexible,
    captionpos=b,   
}

\title{Deal.II 求解泊松方程简化形式}
\author{张志心\\Mixed Class 3210106357}

\begin{document}
\maketitle
\tableofcontents
\newpage
\CTEXsetup[format={\Large\bfseries}]{section}

\begin{abstract}
本项目为本人的数学软件项目作业，复现 deal.II 框架的 examples 中的 step-3 的程序，并分析其结果。deal.II是一个开源的有限元软件库，专门用于求解各种偏微分方程和多物理场问题。step-3 中的程序基于前两个例子，实现了一个拉普拉斯求解器，用于求解柏松方程的简化版本，通过该例，可以初步了解有限元程序的一般结构，包括组装线性系统，求解这个线性方程的方法。本项目还对输出数据进行了可视化展示。
\end{abstract}
\textbf{关键词：} 有限元方法，柏松方程，数据可视化

\section{项目介绍}
\subsection{有限元方法基本设置}
考虑如下简化版本的Poisson方程，其边界值为0，等式右端非0\cite{dealii-step3}：
\begin{equation*}
\begin{aligned}
    -\Delta u &= f ~~~~~ \text{in ~ } \Omega \\
    u &= 0 ~~~~~ \text{on ~ } \partial \Omega
\end{aligned}
\end{equation*}

在 $\Omega = [-1,1]^2$ 这个正方形上求解该方程，令 $f(x)=1$。利用有限元方法求解近似解 $u$，需要首先推导出上述方程的弱形式，通过将方程左乘一个测试函数 $\phi$，并且在域 $\Omega$ 上积分：

\begin{equation*}
    -\int_\Omega \phi \Delta u = \int_\Omega \phi f
\end{equation*}

通过分部积分将上式转化为：

\begin{equation*}
    \int_\Omega \nabla \phi \cdot \nabla u - \int_{\partial \Omega} \phi \mathbf{n} \cdot \nabla u = \int_{\Omega} \phi f.
\end{equation*}

测试函数 $\Phi$  必须满足边界条件（即其来自解的集合的切空间），因此在边界上 $\Phi = 0$，上述方程的弱形式为：

\begin{equation*}
    (\nabla \phi, \nabla u) = (\phi,f),
\end{equation*}

这里我们记 $(a,b) = \int_{\Omega} ab$。

然后，问题就变成了从适当的空间（这里是 $H^1$）中找出一个函数 $u$，其对所有测试函数 $\Phi$ 均成立。一般情况下无法使用计算方法找到这样的函数，而是寻找一个近似的 $u_h(\mathbf{x}) = \sum_j U_j\phi_j(\mathbf{x})$，其中 $U_j$ 是一个待定的未知膨胀系数（自由度），$\phi_j(\mathbf{x})$ 是将使用的有限元形状函数。其中形状函数的定义依赖以下四个内容：

\begin{enumerate}
    \item 用于定义形状函数的网格；
    \item 用于描述在参考单元上使用形状函数的有限元（在\texttt{deal.II} 中，它总是单位线段 $[0,1]$，单位正方形 $[0,1]^2$，或者单位立方体 $[0,1]^3$）；
    \item \texttt{DoFHandler} 对象，以有限元对象提供的参考单元描述为基础，枚举网格上所有自由度；
    \item 从参考单元上的有限元类定义的形状函数到实际单元上的形状函数的映射。
\end{enumerate}

根据以上过程得到一组函数 $\phi_i$，由此定义离散问题的弱形式：找到一个函数 $u_h$（即找到膨胀系数 $U_j$），使得：

\begin{equation*}
    (\nabla\phi_i,\nabla u_h) = (\phi_i, f),~~~i = 0,1,\cdots,N-1.
\end{equation*}

若将 $u_h$ 表示为 $u_h(\mathbf{x}) = \sum_j U_j\phi_j(\mathbf{x})$， 则方程可以写为如下线性形式：

\begin{equation*}
    \begin{aligned}
        (\nabla \phi_i, \nabla u_h) &= \bigg(\nabla \phi_i, \nabla \big[\sum_j U_j \phi_j \big]\bigg) \\
        &= \sum_j \big(\nabla \phi_i, \nabla [U_j \phi_j]\big) \\
        &=\sum_j \big(\nabla \phi_i, \nabla \phi_j\big)U_j
    \end{aligned}
\end{equation*}

问题等价于找到向量 $U$，满足
\begin{equation*}
    \begin{aligned}
        &AU = F\\
        &\text{其中，}A_{ij}=(\nabla \phi_i, \nabla \phi_j) \\
        &~~~~~~~~~F_i = (\phi_i, f)
    \end{aligned}
\end{equation*}

\subsection{$A_{ij},F_i$计算方法}

在 \texttt{deal.II} 中，$A$ 的实例是 \texttt{SparseMatrix} 类型，$U,F$ 的实例是 \texttt{Vector} 类型。在有限元方法中，一般采用\textbf{正交法}计算积分，首先把 $\Omega$ 上的积分分为所有单元上的积分，然后使用每个单元上的一组正交点的加权和来计算。即：

\begin{equation*}
    \begin{aligned}
        A_{ij} &= (\nabla \phi_i, \nabla \phi_j) = \sum_{k\in \mathbb{T}} \int_K \nabla \phi_i \cdot \nabla \phi_j \\
        F_i &= (\phi_i, f) = \sum_{k\in\mathbb{T}}\int_K \phi_i f
    \end{aligned}
\end{equation*}

然后使用正交法计算每个单元上积分的近似值：

\begin{equation*}
    \begin{aligned}
        A_{ij}^K &= \int_K \nabla \phi_i \cdot \nabla \phi_j \approx \sum_q\nabla \phi_i(\mathbf{x}_q^K)\cdot \nabla \phi_j (\mathbf{x}_q^K) w_q^K \\
        F_i^K &= \int_K \phi_i f \approx \sum_q \phi_i (\mathbf{x}_q^K) f(\mathbf{x}_q^K) w_q^K
    \end{aligned}
\end{equation*}

其中 $\mathbb{T} \approx \Omega$，$x_q^K$ 是单元 $K$ 上第 $q$ 个正交点，$w_q^K$ 是其正交权重。

\subsection{求解线性系统}

采用共轭梯度法（Conjugate Gradient method）求解线性系统的解。CG 方法适用于对成正定的矩阵，并且很适合处理大型稀疏矩阵。而且由于有限元方法中只需要求解精确解，所以不需要经过 $N$ 次迭代来计算，并且该方法的复杂度为 $O(N^2)$，比高斯消元方法效率高很多。

在CG方法中，要求解 $Ax = b$ 的近似解 $\widetilde{x}$ 时迭代的终止条件为 $\Vert A(x-\widetilde{x}) \Vert \le\tau$，一般情况下设 $\tau = 10^{-6} \Vert b \Vert$。

\subsection{实现方法}

程序主类如下：

\begin{lstlisting}[language=C++]
class Step3
{
  public:
    Step3 ();
    void run ();
 
  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;
 
    Triangulation<2>     triangulation;
    FE_Q<2>              fe;
    DoFHandler<2>        dof_handler;
 
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
};
\end{lstlisting}

程序需要一个 \texttt{Triangulation} 和一个 \texttt{DoFHandler} 对象，以及一个描述想要使用的各种形状函数的有限元对象。然后还需要一个系统矩阵和右向量以及解向量，以及一个矩阵稀疏的对象。 \texttt{FEValues} 对象只在整个装配过程中需要，因此作为局部变量。

包含的成员函数如下：
\begin{itemize}
    \item \texttt{make\_grid()} ：预处理函数，它设置了 \texttt{triangulation} 对象；
    \item \texttt{setup\_system()}：用于设置问题中的所有数据结构，它的功能是初始化 \texttt{DoFHandler} 对象并且设置线性系统的各个参数。
    \item \texttt{assemble\_system()}：计算 $A_{ij}$ 和 $F_i$ 。
    \item \texttt{solve()}：求解线性系统。
    \item \texttt{output\_results()}：用于处理可视化输出。
\end{itemize}

上述程序均为 public 公有函数，并通过 \texttt{Step3::run()} 来依次调用它们：

\begin{lstlisting}[language=C++]
void Step3::run()
{
  make_grid();
  setup_system();
  assemble_system();
  solve();
  output_results();
}
\end{lstlisting}

程序的主函数为：

\begin{lstlisting}[language=C++]
int main()
{
  deallog.depth_console(2);
 
  Step3 laplace_problem;
  laplace_problem.run();
 
  return 0;
}
\end{lstlisting}

\section{项目结果}

\subsection{问题重述}

本项目在 $[-1,1]^2$ 上求解柏松方程的简化形式：

\begin{equation*}
    \begin{aligned}
        -\Delta u &= 1 ~~~~~ \text{in ~ } [-1,1]^2\\
        u &= 0 ~~~~~ \text{on ~ } \partial [-1,1]^2
    \end{aligned}
\end{equation*}

\subsection{参数设置}

利用上述 Deal.II 框架进行求解。

\begin{itemize}
    \item 在创建网格的时候，对所有单元格进行了 5 次细化，最终网格为 $32\times 32$ 个单元，共 1024 个。
    \item 设置多项式次数为 1（双线性元素），自由度为形函数系数的数量，在线性系统中，每个顶点对应一个形函数，分别只有一个系数，自由度为1。在 $32\times 32$ 网格上一共有 $33 \times 33$ 共 1089 个顶点，因此自由度 (degrees of freedom) 为 1089。
    \item 设置 CG 算法迭代次数上限为 1000，迭代终止条件为 $\tau = 10^{-6} \Vert b \Vert$，其中 $b$ 为线性系统右端向量。
\end{itemize}

\subsection{输出结果}

程序运行输出如下：

\begin{verbatim}
    Number of active cells: 1024
    Number of degrees of freedom: 1089
    DEAL:cg::Starting value 0.121094
    DEAL:cg::Convergence step 36 value 7.07761e-08
\end{verbatim}

说明 CG 求解器进行了 36 次迭代之后终止，第三行行第四行的值分别为开始和终止时的残差的范数。程序运行同时还生成了 \texttt{solution.vtk} 文件，是解函数的三维可视化文件，使用 visIt 软件处理之后得到 \texttt{*.bmp} 图片如下所示：

\begin{figure}[H]
  \centering
  \begin{tikzpicture}
    \node[inner sep=0pt] (image) at (0,0) {\includegraphics[width=0.75\textwidth, bb=0 0 36cm 33cm]{../result/protected/solution10001.bmp}};
    \draw [black] (image.south west) rectangle (image.north east); 
  \end{tikzpicture}
  \label{fig:solution.vtk}
  \caption{solution.vtk 可视化}
\end{figure}

泊松方程等价地描述了一个膜在收到外部垂直力之后变形的形态，恒定的压强下，膜会自然向上鼓起，就和上图所显示的一样。

\bibliographystyle{plain}
\bibliography{ref}


\end{document}